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In situations involving large potentials or surface charges, the Poisson Boltzman(PB) equation 
has shortcomings because it neglects ion-ion interactions and steric effects. This has been widely 
recognized by the electrochemistry community, leading to the development of various alternative 
models resulting in different sets " modified PB equations" , which have had at least qualitative suc- 
cess in predicting equilibrium ion distributions. On the other hand, the literature is scarce in terms 
of descriptions of concentration dynamics in these regimes. Here, adapting strategies developed 
to modify the PB equation, we propose a simple modification of the widely used Poisson-Nernst- 
Planck (PNP) equations for ionic transport, which at least qualitatively accounts for steric effects. 
We analyze numerical solutions of these MPNP equations on the model problem of the charging of a 
simple electrolyte cell, and compare the outcome to that of the standard PNP equations. Finally, we 
repeat the asymptotic analysis of Bazant, Thornton, and Ajdari 1] for this new system of equations 
to further document the interest and limits of validity of the simpler equivalent electrical circuit 
models introduced in Part I for such problems. 



I. INTRODUCTION 

In Part I of this series, we focused on steric effects 
of finite ion size on the charging dynamics of a quasi- 
equilibrium electrical double layer, motivated by the 
breakdown of the classical Gouy-Chapman model 0] at 
large voltages (up to several Volts and 3> kT/e = 25 
mV), e.g. which are commonly applied in AC electro- 
osmosis [1, 0, H, H, 01 ■ We introduced two simple mod- 
ifications of the Boltzmann equilibrium distribution of 
ions to incorporate steric constraints. Both new models 
predicted similar dramatic consequences of steric effects 
at large voltages, such as greatly reduced diffuse-layer 
capacitance and neutral salt uptake from the bulk com- 
pared to the classical theory. The crucial effect of the 
finite ion size is to prevent the unphysical crowding of 
point-like ions near the surface at large voltages by form- 
ing a condensed layer of ions at the close-packing limit 
(likely to include at least a solvation shell around each 
ion). 

The idea that the electric double layer acts like a ca- 
pacitor is part of a bigger picture and suggests that the 
dynamics can be described in terms of equivalent cir- 
cuits [1, where the double layer remains in quasi- 
equilibrium with the neutral bulk. This classical ap- 
proximation has been discussed and validated in the 
thin double layer limit by asymptotic analysis of the 
Poisson-Nernst-Planck (PNP) equations [H. The PNP 
equations provide the standard description of the linear- 
response dynamics of electrolytes perturbed from equilib- 
rium, based on the same assumption of a dilute solution 
of point-like ions interacting through a mean field which 
underlies the PB equation for equilibrium [3, fiot . 

Here we try to account for the effect of steric con- 
straints on the dynamics, by first establishing modified 



PNP (MPNP) equations using linear response theory 
and modified electrochemical potentials. We apply the 
MPNP equations to describe the charging of a parallel 
plate electrolyte cell in response to a suddenly applied 
voltage and comment on the differences with usual PNP 
dynamics. The results are in line with the work in Part I 
since the double layer behaves like a capacitor; however 
its capacitance is reduced by steric effects, and neutral 
salt uptake is decreased as well. 

Finally, following the analysis of Bazant, Thornton, 
and Ajdari we demonstrate that the considerations of 
Part I can rigorously be supported by asymptotic analysis 
on the MPNP equations. This helps us understand the 
limits of the electric double layer capacitor models and 
define rigorously what is meant by the thin double layer 
limit. 



A. Electrolyte Dynamics in Dilute Solution Theory 

As we have mentioned in Part I, the dilute solution 
theory has been the default model for ion transport for 
the most part of the twentieth century. According to this 
theory, it is acceptable to neglect interactions between 
individual ions. As a result, the electrochemical potential 
takes the form 

= fcThiQ + z,e7/' (I) 

where z^e is the charge, Ci the concentration, and ■0 the 
electrostatic potential-usually governed by the Poisson's 
equation. In order to derive an equation for the elec- 
trolyte dynamics, we need to combine the above equation 
with the flux formula (here with the standard assumption 
that interactions between different species are negligible. 
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see e.g. 



(2) 



(where 6j is the mobihty of the species i) and the general 
conservation law, 



dcj 
dt 

The resulting equations, 
da 



dt 



(3) 



(4) 



are called the Nernst-Planck (NP) equations. Here, the 
Einstein's relation, Di = bi/kT, relates the ions mobil- 
ity bi to its difFusivity Di. The system is closed by the 
Poisson's equation. 



(5) 



where e denotes the permittivity of the electrolyte. The 
name Poisson-Nernst-Planck(PNP) equations is coined 
for the system defined by equations @ and The 
standard PNP system presented above has been used to 
model selectivity and ionic flux in biological ion chan- 
nels [H, [l^ [H, {TtL [isi. AC pumping of liquids over 

electrode arrays rams, 

l2ll |22| |. induced-charge 
electro-osmotic flows around metallic colloids '2§, and 
microstructures [1, 113,125], dielectrophoresis 26, 27, 28] 
and induced-charge electrophoresis [28|, |29|, (sfl, Ej] of po- 
larizable particles in electrolytes. 

As explained thoroughly in Part I, however, the dilute 
solution theory, including the PNP equations, has limited 
applicability: its predictions easily violate its basic as- 
sumption (i.e. being dilute) near surfaces of high poten- 
tial. In fact, this happens more often than not due to the 
exponential (Boltzmann-type) dependence of counter-ion 
concentration on electrostatic potential. The steric limit, 
that is Cmax = l/a"^, a being the typical spacing between 
densely packed ions, is reached at the critical potential 



-— ln(coa3) ^—\n{^ 

ze ze Co 



(6) 



where cq is the bulk electrolyte concentration of either 
of the species. Due to the logarithmic dependence in its 
formulation, the critical voltage ^'c is no more than a 
few times the thermal voltage ipT = kT / ze and therefore 
easily reached in many applications such as the induced 
charge electroosmosis 23, 24]. This has motivated re- 
searchers to modify the standard equations and improve 
the dilute solution theory (see next subsection). 



B. Beyond Dilute Solution Theory 

Statistical mechanics have proved to be an indispens- 
able tool for analyzing and improving the dilute solution 



theory. In particular, a statistical model with the desired 
level of detail can be set up, and after the corresponding 
free energy F is calculated, the chemical potentials can 
be obtained from the formula 



Ml 



(7) 



where ci is the concentration of the i-ih species. Then 
differential equations governing the electrolytes can eas- 
ily be derived from this chemical potential as outlined in 
the previous subsection. There have been many attempts 
to calculate the free energy of the electrolyte more accu- 
rately to improve the dilute solution theory (see below for 
references). In some of these attempts, the calculation of 
the free energy was replaced by an equivalent consider- 
ation of mean electrostat ic p otential and mean charge 
density [33,|3l,[3l,|3l,[3l,[33- Using those mean quanti- 
ties, various correlation functions as well as new and more 
accurate PB (i.e. MPB) equations have been proposed. 
However, one should keep in mind that the correspond- 
ing free energies can still be calculated for these models, 
too. 

Perhaps the first examination of the limits of the dilute 
solution theory by statistical mechanical considerations 
was by Kirkwood [38] in 1934. After a detailed analysis 
of the approximations of the dilute solution theory, Kirk- 
wood concluded that those approximations consisted of 
the neglect of an exclusion- volume term and a fluctuation 
term. Furthermore, he gave estimates of those terms and 
argued that they are indeed negligible in the bulk elec- 
trolyte. In recent years, there have been many attempt 
[3^, [33I , [3^ . [ssl . [sgI . Wf\ originating from the liquid-state 
theory to calculate those neglected terms more explicitly, 
which resulted in a variety of MPB equations (includ- 
ing the MPB1,...,MPB5,:. hierarchy of Outhwaite and 
Bhuiyan H). 

Another general approach to the statistical mechan- 
ics of electrolytes is based on density functional theory 
(DFT). Using this formalism, Rosenfeld [3l|40,|4l| sys- 
tematically derived elaborate free energy functionals for 
neutral and charged hard-sphere liquids starting from ba- 
sic geometric considerations. Gillespie et al. [i^, [11, [il] 
calculated chemical potentials from Rosenfeld's free en- 
ergy functionals, and used them along with the formula 
to calculate the steady flux in an ion channel, as well 
as to investigate the equilibrium structure of the double 
layer. With this theoretical framework, Roth and Gille- 
spie [4^ were able to explain size selectivity of biological 
ion channels. 

Perhaps one reason why neither the MPB hierarchy of 
Outhwaite and Bhuiyan (33| . nor the free energy func- 
tionals of Rosenfeld [4l[ , has gained widespread use and 
recognition is their intrinsic complexity which limits their 
simple application to specifc problems. For example, the 
hard sphere component of Rosenfeld's free energy density 
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is given by 

^HS = -noln(l - n^) 

nin2 - nvinv2 nj f. _ "^^2^ 

1 - ns 247r(l - 713)2 J 

where each of no, ^i, 722, 723, nyi, and 77^2 are functions 
defined in terms of (in general 3D) integrals. So the im- 
provements over the initial PB equations are made at 
the cost of losing the possibility of analytical progress, 
and nontrivial numerical work is required even for sim- 
ple problems in one dimension. 

A considerably simpler approach, which mainly focuses 
on the contribution of size effects to the free energy, is 
based on mean-field theory together with the lattice-gas 
approximation in statistical mechanics. Following the 
early work of Eigcn and Wicke [H, 113, HH, Iglic and 
Kral-IgliciiJsaM 112 and Borukhov, Andelman, and 
Orland [53l. I54L l55l| were able to come up with a sim- 
ple statistical mechanical treatment which captures ba- 
sic size effects. A free energy functional is derived by 
mean-field approximations of the entropy of equal-sized 
ions and solvent molecules. This free energy is then min- 
imized to obtain the equilibrium average concentration 
fields, and the corresponding modified PB (MPB) equa- 
tion. While more complex than the original PB equa- 
tion, it is still rather compact and simple, and definitely 
amenable to further analytical use. 

C. Scope of the present work 

All of the above authors, as well as others we have cited 
in Part I, focus on the equilibrium or steady state prop- 
erties of the electrolytes. In fact, we are not aware of any 
attempt to go beyond the dilute solution theory (PNP 
equations) in analyzing the dynamics of electrolytes in 
response to time-dependent perturbations, such as AC 
voltages. Here in the second part of this series, our aim 
is to improve the (time-dependent) PNP equations of the 
dilute solution theory by incorporating the steric effects 
in a simple way with the goal of identifying new generic 
features. Our focus is electrolyte systems that contain 
highly charged surfaces, such as an electrode applying a 
large voltage V ^ tpT = kT/ze , where the steric effects 
become important quantitatively as well as qualitatively. 
As in Part I, here we a gain adopt the mean-field ap- 
proach of Borukhov et al.[53| to the size effects, because 
of two main reasons: First, and foremost, it is preferable 
to start with simple formulations that capture the essen- 
tial physics while remaining analytically tractable, as we 
are mainly interested in new qualitative phenomena. Sec- 
ond, it is not clear to us how well the liquid-state theories 
would perform at large, time-dependent voltages, since 
they are (at least in some cases) based on perturbative 
methods around an equilibrium reference state -which 
may not even exist, say, in presence of an AC electric 
field. 



A quick outline of the paper is as follows: In sec- 
tion II, we derive the modified Poisson Nernst-Planck 
(MPNP) equations, using the free energy obtained by 
[i^llll to derive the modified PB equations (MPB). In- 
stead of minimizing the free energy F we first compute 
the chemical potentials from the equation ([7]). Then we 
describe the electrolyte dynamics by linear response re- 
lations described by equations ((H), ([3]). As a result, we 
end up with the promised MPNP equations which include 
steric corrections to the standard Nernst Planck equation 
that become increasingly important as the concentration 
field gets large. We continue in section III by setting up 
and investigating the numerical solutions of our modified 
PNP equations for the problem of parallel plate block- 
ing electrodes. In section IV, we follow the same lines 
as in and establish our earlier conclusions (includ- 
ing the electric circuit picture) about electrical double 
layers in Part I by a rigorous asymptotic analysis. We 
also calculate higher order corrections to the thin dou- 
ble layer limit, and check a posteriori the validity of the 
leading order approximation. Finally, we close in section 
V by some comments and possible future directions for 
research. 



II. DERIVATION OF MODIFIED PNP 
EQUATIONS 

In this section, we derive the modified PNP equations 
as outlined in the introduction. For simplicity, let us re- 
strict ourselves to the symmetric z : z electrolyte case. 
We also assume that the permittivity e is constant in the 
electrolyte. In the mean-field approximation, the total 
free energy, F — U — TS, can be written in terms of 
the local electrostatic potential ip and the ion concentra- 
tions c± . Following [53] , we write the electrostatic energy 
contribution U as 

U = J dr (^-^iVi^f + zec+^p- zec-tf^ (8) 

The first term is the self-energy of the electric field for a 
given potential applied to the boundaries (which acts as 
a constraint on the acceptable potential fields), the next 
two terms are the electrostatic energies of the ions. The 
entropic contribution (the steric effects) can be modeled 
as [5jj 

- TS* / dr{c+a^ In (c+a^) + c^a^ In (c^a^) 

+ (1 - c+a^ - c-a^) In (l - c+a^ - c-a^)} (9) 

where we have assumed for simplicity that both types of 
ions and solvent molecules have the same size a. The first 
two terms are the entropies of the positive and negative 
ions, whereas the last term is the entropy of the solvent 
molecules. It is this last term, which penalizes large ionic 
concentrations. 

Requiring that the functional derivatives of this free 
energy F with respect to 7/; and c± be respectively zero 



and constant chemical potentials /z± (the Lagrange mul- 
tipliers for the conserved number of particles of each 
kind), Borukhov et al.[5^ obtain the modified PB equa- 
tion 



zeco 



2 sinh 



[kT ) 



1 



2j/sinh2 (ir) 



(10) 



Here we go one step further and derive MPNP by calcu- 
lating the chemical potentials fi± from yielding 



= C (11) 
— ±zeip + kT{\nc± — In (l — c+a^ — c_a'^)) 



and as a reasonable form for the dynamics we postulate 



dc± 
'd7 



V■{b±c±V^^±) 



(12) 



This form which is classical for close to equilibrium trans- 
port is already an approximation, as it neglects cross- 
terms in the mobility matrix (i.e. that the gradients in 
can induce a current of positive ions). Further, we 
now assume that the mobilities for each type of ions are 
the same, and equal to 6 = 5+ = &_ which is consistent 
with the assumption that they have the same effective 
size a. A final approximation is that we take this value b 
to be constant, and in particular insensitive to the crowd- 
ing that can occur in the electric double layers. All these 
approximations will be further discussed and challenged 
in subsequent work, but we proceed for the time being 
with the present simpler version, which is a first attempt 
to incorporate steric effects in the dynamics. As the elec- 
tric fields adjusts almost instantaneously to minimize the 
electrostatic energy, we get Poisson equation as 



+ ze (c+ — c_) =0 



(13) 



The equations (I12p yield modified Nernst-Planck equa- 
tions: 



dc± 
'd7 



D 



DV^c±±- zeV • (c± Vi/-) + a-'DV ■ ( 

kbT 



(14) 
c±V (c+ -I- c_) 



1 — c+a'^ 



) 



where D = kTb is the common diffusion coefficient. 
As mentioned before, the extra last term is a correc- 
tion due to the finite size effects. Considered together 
the above set of equations are modified Poisson-Nernst- 
Planck equations (MPNP), which is our simplest pro- 
posal for a dynamic description that incorporate steric 
effects. As in the standard case, this set of equations are 
completed by appropriate boundary conditions. In par- 
ticular, we consider in this paper that no reactions take 
place at the surface (blocking electrodes) so that no-flux 
boundary conditions 



DWc± ± bzec±Vtp 



a?Dc±\/ (c+ + C-) 
1 — c+a-^ — C-O^ 







(15) 
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FIG. 1; Sketch of the model problem (from Ref. Ql). A 
voltage 2V is suddenly applied to a dilute, symmetric, binary 
electrolyte between parallel-plate, blocking electrodes sepa- 
rated by 2L. 



hold for the ions. We write the boundary condition for 
the potential by accounting as before for the possible 
presence of a thin insulating layer of fixed capacitance 
Cs- This leads to a mixed boundary condition 

dip 

'(/)(n = 0) = Velectrode + AsT^(n = 0) (16) 

on 

where ipeiectrode is the applied potential at the electrode 
which is reduced by the insulating layer to il'n=o a-t the 
surface of the electrolyte (where MPNP starts to be ap- 
plied). Xs = e/Cs is a measure of the thickness of this 
layer. Here, n is the normal direction to the surface point- 
ing into the electrolyte. 



III. MODEL PROBLEM FOR THE ANALYSIS 
OF THE DYNAMICS 

In order to gain some insight into the ramifications of 
the extra term we introduced into PNP, we now turn 
back to the basic model problem discussed in reference 
Namely, as shown in Fig. [U we consider the ef- 
fectively one-dimensional problem of an electrolyte cell 
bounded by two parallel walls (at x = ±L), filled with 
a z:z electrolyte, at concentration Cq, and across which 
a step voltage of amplitude (21^) is suddenly applied at 
t = 0. We further assume that no Faradaic reactions are 
induced at the electrodes surface. We formulate MPNP 
in this setting, and then compare some numerical solu- 
tions of MPNP to those of PNP. In the next section we 
will focus on the case where the double layer thickness 
Xd = (2e^co/£fcBr)~^/^ is much smaller than the (half) 
cell thickness L. 

First off, we note that, in this simple geometry, the 
gradients are replaced by ^ and the derivative with re- 
spect to surface normal is replaced by ^ at a; = — L, 
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and by - at x = L. Following fi\, we cast the MPNP 
equations in dimensionless form using L as the reference 
length scale, and Tc — XdL/D as the reference time scale, 
thus time and space are represented by i = tD/X]jL and 
X = x/L. The problem is better formulated through the 
reduced variables c = (c+ + c_ ) for the local salt con- 
centration, p — 2^ (c-|_ — c_) for the local charge density, 
and (p = zetp/kT for the electrostatic potential. The so- 
lution is determined by three dimensionless parameters: 
V = zeV/ kT, the ratio of the applied voltage to thermal 
voltage, e = Xd/L, the ratio of the Debye length to the 
system size, 6 — Xs/ Xd introduced in section II which 
measures the surface capacitance, and finally a new pa- 
rameter quantifying the role of steric effects v = 2a^co, 
the effective volume fraction of the ions at no applied 
voltage. After dropping the tildes from the variable x, 
the dimensionless equations take the form 



dc 






50 


dt 


dx ' 


\dx 




dp 


= e— 1 


fdp 


d<i> 


dt 


dx ' 


\dx 


dx 








d^cj) 








dx^ 



dc 



1 — vcdx 
vp dc 
I — vc dx 



(17) 



with completely blocking boundary conditions at x = ±1, 



5^ 
^ dx 



dc 
dx 

dp ^ d(j) 

dx dx 



dc 



1 — vcdx 
p dc 



1 — vcdx 
in addition to which reads 



' ± (Se— — — ±v 
dx 



= 




(18) 



(19) 



at a; = ±1, where again 5 measures the effective thickness 
of the surface insulating layer. Because it is impossible to 
satisfy all the boundary conditions when e = 0, the limit 
of vanishing screening length, e ^ 0, is singular. The 
total diffuse charge near the cathode, scaled by 2zecoL, 
is 



lengthy, we focus on a few situations for which we pro- 
vide plots meant to illustrate the differences brought in 
by accounting for steric effects. We therefore also plot the 
outcome of classical PNP in the same situations (which 
correspond to = in the equations). Of course a sys- 
tematic difference is that with the MPNP neither the 
concentration c nor the charge density p ever overcome 
the steric limit I /v. 

For sake of readability of the figure, we start in Fig[2] 
with untypically large values for both the e = 0.1 and 
I' = 0.25. The potential is ten times larger than the 
thermal voltage u = 10. The map of the concentration 
and charge density are given for different instants after 
the application of the potential drop. Build-up of the 
double layers, and the consequent depletion of salt in 
the bulk are visible. The MPNP solution stays bounded 
by l/v as promised, whereas the PNP solution blows up 
exponentially. A consequent observation is that salt de- 
pletion in the bulk is weaker with the MPNP, whereas 
with the classical PNP bulk concentrations drop to small 
values even for this moderate potential (^ 0.25 V in di- 
mensional units). Of course, this is also a consequence of 
the large e. 

For a simulation with more realistic values, we have 
taken v — AO (corresponding to 1 in dimen- 
sional units), and e — — 0.01. The correspond- 
ing solutions are plotted at nondimensional times t = 
0, 0.5, 1, 2, 4, 6, 8, 10, ...,50 in FigEl Figure ^Si) shows 
the bulk concentration dynamics, whereas Figure[3](b) fo- 
cuses on the double layer near the boundary at x = — 1. 
The MPNP solution again stays bounded by 1/1/ = 100. 
The PNP solution is not plotted as it blows up in the 
double layer to about cosh(40) « 10^^ times the bulk 
value (itself consequently very small), requiring subtle 
numerical methods. With the MPNP, the charge build 
up of the double layer first proceeds as it would with the 
PNP until concentrations close to the threshold l/v are 
reached. Thereafter, charging proceeds by growth of the 
double layer thickness at almost constant density. 



IV. ASYMPTOTIC ANALYSIS 



(20) 



The dimensionless Faradaic current, scaled to 2zecQD/ L, 
is 



dp ^ d(f) ^ p dc 
dx dx 1 — vc dx 



(21) 



We have numerically solved these equations for various 
values of the dimensionless parameters. As a complete 
description of this large parameter space would be very 



In this section, we adapt to the MPNP equations in- 
troduced here the asymptotic analysis presented in [l[ 
for PNP equations in the limit of thin double layers. We 
show rigourously that the key properties of the charg- 
ing dynamics remain the same, namely, at leading order 
the boundary layer acts like a capacitor with a total sur- 
face charge density q{t) that changes in response to the 
Ohmic current density j(t) from the bulk. However, the 
expressions for the capacitance of the double layer, and 
the diffusive flux from the bulk into the double layer do 
change, as we already discussed in Part I. 
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V = 0.01, e = 0.01, ■u = 40 




-0.5 



0.5 




(b) 

FIG. 2: The numerical solutions to the PNP and MPNP sys- 
tems. The dimensionless charging voltage is w = 10, which 
corresponds to approximately 0.25 Volts at room tempera- 
ture. The very large values of e — 0.1 and v — 0.25 are 
chosen deliberately for illustration. The dimensionless bulk 
concentration field c is shown in (a), and the dimensionless 
charge density p in (b) 




V = 0.01, e = 0.01, D = 40 



100 



t = 0, 0.5, 1, 
2, 4, 6, 8, 10, 
12, 50. 




(b) 



-0.998 -0.996 -0.994 -0.992 -0.99 

X 



FIG. 3; The numerical solution to the MPNP system at the 
dimensionless times (i.e. time scaled to the charging time) 
t — 0, 0.5, 1, 2, 4, 6, 8, 10, 50. The dimensionless charging 
voltage is w = 40, which corresponds to approximately 1 Volt 
at room temperature. The parameters t = v = 0.01 are still 
large, but comparatively realistic. The dimensionless bulk 
concentration field c is shown (a) globally for the whole re- 
gion (b) zoomed in at the double layer. 



therefore, we consider only — 1 < a; < 0. 

The outer solution, in the bulk, is denoted by a bar 
accent, and the expansion takes the form 



A. Inner and Outer Expansions 

c = c{x, i) = Co + eci + ... 

Adopting the notation in we seek regular asymp- 
totic expansions in e. We observe that the system has jt is easy to check that cq = 1, po = and (^o = as 
the following symmetries about the origin ^ direct consequence of our choice for the time scale as 

Tc. As to the inner solutions (i.e. a; w — 1), we remove the 
c(— a;, i) = c(x, t), p(— a;, i) = — p(a;, i), 4'{—x,t) — —(f>{x,t) singularity of equation ([T7]) by introducing ^ = (1 -I- a;)/e. 
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Then, we obtain 



dc d f dc dcj) vc dc 
dp d f dp 90 vp dc 



(22) 
(23) 
(24) 



dt d^ \di d^ l~iycd^ 

for which, we can seek regular asymptotic expansions 
c = c(^, t) = Co + eci + ... 



Matching inner and outer solutions in space involves the 
usual van Dyke conditions, e.g. 

lim c(^, t) = lim c{x, t) 

which implies co(oo,i) = co(— ci(oo,i) — ci(— 1,<), 
etc. As seen from equations ([^ - (1^ . there are no terms 
with time derivatives at leading order. This quasiequilib- 
rium occurs, because the charging time Tc is much larger 
than the Debye time td = Xjj/D, which is the character- 
istic time scale for the local dynamics in the boundary 
layer. Consequently, the leading order solution is the 
" equilibrium" 



Co 



cosh $0 



Pa 



— sinh $0 



l + i/(^cosh$o- Ij l + z/(cosh$o- 1) 

where the excess voltage relative to the bulk 

$ = 0(e, t) ~ ^(-1, t) - l-o + el>i + ... 
satisfies the modified PB equation at leading order 



d^ 

de 



sinh $0 



l + u [ cosh $0-1 



(25) 



with $0(00, t) — 0, and $0(0, i) = C(*): the dimensionless 
zeta potential, which varies as the diffuse layer charges. 
Applying matching to the electric field, we obtain 

|(oo,t)^.|(-l,t)^ f (oo,.)=0 (26) 
and therefore an integration of ((25)) yields 



^ = -sign(|.o)y^ In ( f + I. ( cosh $0 - 1 ) ) (27) 



B. Time-dependent matching 

So far we have found that both the bulk and the bound- 
ary layers are in quasiequilibrium, which apparently con- 
tradicts the dynamic nature of the charging process. This 



is reconciled by noting once more that the boundary con- 
dition for the inner solution involves the quantity C,{t), 
which varies in response to the diffusive flux from the 
bulk. Motivated by the physics, we consider the total 
diffuse charge, which has the scaling q{t) ~ cq[t), where 



j'OO 

Jo 



qo + eqi + f^q2 



(28) 



Taking a time derivative, and using (j23p together with 
the no flux boundary condition (jlSp . we obtain 



i^p 



dc 



lim ( ^ 

E-f-i V dx 



1-i^cdC 
vp dc 



dx 1 — vc dx 



where we applied matching to the flux densities. Substi- 
tuting the regular expansions of inner and outer solutions 
yield a hiararchy of matching conditions. At leading or- 
der, we have 



dt 



(29) 



which, being a balance of 0(1) quantities, is reassuring 
us that we have chosen the correct time-scale. This is 
the key equation which tells us that at leading order, the 
double layer behaves like a capacitor, whose total surface 
charge density g, changes in response to the transient 
Faradaic current density, j(t), from the bulk. 



C. Leading Order Dynamics 



Using equations ([241) , (US]) and (|27l), the integral in ([28 
can be performed at leading order to yield 



qo{t) - -sign(Co)y- ln[l + z.(coshCo - 1)] (30) 
Then the Stern boundary condition, ea. p^ . yields 



Co + feign(Co; 



■ln[l + i/(coshCo - 1) 



(31) 



where *(t) = -v - (j){-l,t) ^ + e#i + ... is the 
total voltage across the compact and the diffuse layers. 
Equation ([3T|) results in higher Co than its classical PNP 
(or PB) counterpart 



Co + 2,5sinh(Co/2) = *o 



(32) 



because the left hand side of (|3T|) is always smaller than 
the left hand side of ([5^ for any given Co- In both for- 
mulas, the first term (i.e. Co) is the voltage drop over 
the diffuse layer whereas the second term (i.e. —Sqo) is 
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the voltage drop over the Stern layer. For high applied 
voltages, the standard model assigns exponentially big- 
ger proportions of that voltage to the Stern layer, whereas 
the modified theory predicts a balanced distribution. For 
more details, see Part I section V. 

Substituting into the matching condition (|29|1. we ob- 
tain an ordinary initial value problem 

-Coi^o)^^^o + v, MO)=v (33) 



where 



dqn 

- — 



(34) 



l + i/(coshCo-l) / 2 



I sinh Co I 



^ln(l -Hi/(coshCo - 1)) +<5 



is the differential capacitance for the double layer as a 
function of its voltage. This has a completely different 
behaviour than its counterpart in namely 



Co 



1 



sech(Co/2) + (5 



(35) 



especially at higher 



As Co 



oo, our differential ca- 
pacitance Cq ~ {J ^Co + S)^^ 0; in contrast to Cq^°, 
the differential capacitance given by (f35|) , which tends to 
5""'^ (see Part I for more details). Thus, at high zeta po- 
tentials, steric effects decrease the capacitance, possibly 
down to zero. Physically, this is because of the increasing 
double layer thickness as a result of excessive pile up of 
the ions coming from the bulk into the double layer (see 
Part 1). 

Equation is separable, and its solution can be ex- 
pressed in the form 



^o=Joit)-v = -F-Ht) 



F{z) 



Co{u) 
V + u 



du 



(36) 



(37) 



Steric effects reduce the capacitance, and therefore F{z), 
which by formula (|37p implies a faster relaxation process 
for the voltage difference ^Pq. In other words, the RC 
relaxation time is shortened (see Part I). 



D. Neutral Salt Adsorption by the Double Layer 

A natural consequence of reduced capacitance at the 
electrodes is the reduction of the amplitude of the dif- 
fusion layer in the bulk were the adsorbed salt in the 
double layer is extracted from. We now revisit some of 
the ideas in and recalculate the neutral salt adsorp- 
tion by the double layer. The excess ion concentration. 



c — Co, acquired by the double layers is accounted for by 
the diffusion from the bulk at 0(e) or higher, as diffusion 
is absent at leading order. Following [l|], we introduce 
the variable 'w{t) = ew{i), akin to q{t), which represents 
the excess amount of salt in the double layer: 

w{t)= / [5{^,t)-co{-l,t)]d^^Wo{t)+^Mt) + - 
Jo 



Taking a time derivative, we find 



dl 



—j—[t) = hm - — -f- p— + — ^ 



' dc 

lim — h p-^ 



dc 



(38) 



-i\dx dx I — vc dx ^ 
Since c = 1 + eci + ... , at leading order c ^ 1 and 

f ~ e^, therefore ^ yields 

^(0-^f^(t) = g(-M) (39) 
which involves the new time variable 



t = 



-t = 



(1 - J.) (1 - ,y) T, 



e T T 



scaled to bulk diffusion time, Ts = (1 — v) L^/Z), which 
is slightly different from the time scale given in mi. The 
salt uptake w[t) can be expressed in terms of an integral 



w{t) 



cosh <i> 



1 + (^cosh $ - 1 
cosh $ - 1 




\ + v (^cosh $ - 1 



2 In ( 1 + ( cosh $ - 1 
(40) 



with no obvious further simplification. 

We now proceed to calculate the depletion of the bulk 
concentration during the double-layer charging. Note 
that the new time scale introduced by p9p is the time 
scale for the first order diffusive dynamics in bulk 



dc\ (1 — v) dci d'^ci 
'dt ^ 'i 



dt 



dx^ 



(41) 



As the source is defined by (|39p in terms of gradients, 
an appropriate Green's function can be obtained by tak- 
ing Laplace transforms and using method of images as in 
[ij, which leads to 



ci(x,t) = - ^ dt'G{x,t- t')^-[wo{t'/e)] (42) 



where 



^ oo 



-{x-2m+lY /it 



(43) 
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In the limit e — > 0, the initial charging process at the 
time scale = O (e) is almost instantaneous, which is 
followed by the slow relaxation of the bulk diffusion lay- 
ers. This limit corresponds to approximating the source 
terms in the integral in p2)) to exist only in a small O (e) 
neighborhood of zero, or more explicitly 



limci(x,t) ~ -G{x,t) '^^'m^'^" (7)] 

= -wq (00) G{x, t) (44) 

with Wq (00) = 

cosh$-l {l-v)d^ 
Jo l + i/(^cosh$ ^ ' 

where 



ij ln(l + j/(coshl>- 1)) 

(45) 



/(C) = C + sign(C)5^/ - In {l + u (coshC - 1)). 



(46) 



As expected from the underlying physics, and already 
explained in Part I and illustrated in the plots of sec- 
tion III above, the formula predicts smaller values 
for the the depth of the bulk diffusion than its classical 
PNP {v — 0) counterpart. This difference is more pro- 
nounced at higher v. Equation (|44p is a simple approxi- 
mation that describes two diffusion layers created at the 
electodes slowly invading the entire cell. At first, they 
have simple Gaussian profiles 



c{x, t) ^ 1 



ewo (00) 



(47) 



for e << f << 1. The two diffusion layers eventually col- 
lide, and the concentration slowly approaches a reduced 
constant value, 



c{x, t) ^ 1 — ewQ (00) 



(48) 



for t >> 1, as we expect from the steady-state excess 
concentration from the double layers. 

Of course one expects that at large enough applied 
voltage, the above approximation breaks down, as the 
decrease in bulk concentration becomes significant and 
thus modifies the value of w in the double layer. 



E. Validity of the Weakly Nonlinear 
Approximation 

The first order solution consisting of the variables in- 
dexed by "0" is often referred to as the weakly nonlin- 
ear approximation, whose main feature is that the bulk 
concentrations are constant, namely c ~ cq — 1 and 
p « Pq = 0. The system therefore is characterized only by 
the surface charge q ~ qo, or the double layer potential 



difference 4* « ^Pq, which is governed by the nonlinear 
ODE in equation ([55]) . This corresponds to modeling 
the problem by an equivalent circuit model with vari- 
able capacitance for the double layer, and constant bulk 
electrolyte resistance. 

In order to understand when the weakly nonlinear ap- 
proximation holds, we can compare the size of the next 
order approximation to the leading term. Although not a 
rigorous proof, one may argue that if |eci | is much smaller 
than Co = 1, then leading order term is a good approx- 
imation to the full solution. We will get help from the 
approximations (|Tfl) and ([^5]) to see if this is the case. 

Seen in the light of equation (|48|) . the assumption 
that the first correction is much smaller than the leading 
term requires that Vb = ewQ (00) <C 1, in other words 



cosh <I> — 1 



« 1, 



^^-^n l + ^'Ccosh*-!) ^ll„(i+^(cosh<i>-l)) 

where / is given in (j46[) . After a series of approxima- 
tions, including (1 — i/) « 1, (5 = 0(1), and ( ^ l(or 



V ^ 1), this becomes ^ In (1 + i/coshu) <C 1. If is 

not too close to zero(i.e. z/coshu ^ 1), this is the same 
as le^vjv <C 1. Putting the units back, we obtain 



A 



zeV 



-r2 3 «1 (49) 

For a typical experiment with = 10 nm, L = O.lmm, 
Co = ImM, a = 5A, and at room temperature, this con- 
dition becomes V <^ 188 Volts. 

However, the weakly nonlinear dynamics breaks down 
at somewhat smaller voltages, because the neutral salt 
adsorption causes a temporary, local depletion of bulk 
concentration exceeding that of the final steady state. In 
our model problem, the maximum change in bulk concen- 
tration occurs just outside the diffuse layers at x = ±1, 
just after the initial charging process finishes at the same 
scale t = 0(1) or t — 0(e). Letting t = e, and x — ±1 
in equation (|47p . we obtain the first two terms in the 
asymptotic expansion as 



c(±l,e) = 1 - Y (00) 

At that time, the double layers have almost been fully 
charged , however the bulk diffusion has only had time 
to reach a region of length O (y/e) .So the concentration 

is depleted locally by O =0 {^/e) , which is much 

more than the uniform 0(e) depletion remaining after 
complete bulk diffusion. 

Therefore, in order for the time-dependent correction 
term to be uniformly smaller than the leading term, we 
need 



Wo (00) <€. 1 



(50) 



TT 

By the same approximations as in p9|) , this yields <C 
1, or with the units 



Xn zeV 



y/n La^co kT 



< 1 



(51) 
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The corresponding threshold voltage is smaller than the 
former by a factor of roughly ^ L/Xj:, — e~^. For the 
same set of parameter as above Ad = 10 nm, L = 0.1mm, 
Co = ImM, a = bA, and at room temperature, this con- 
dition gives V <^ 0.033 Volts. Thus according to this cri- 
terion, the weakly nonlinear approximation easily breaks 
down, and one may need to consult to the full MPNP 
system for an understanding of the electrolyte dynamics. 

A more accurate understanding into the condition ([50)1 
is gained by the numerical study of the function w. The 
curves on which y/ew = 1 (i.e. e — w^"^) are plotted in 
FigHlfor various values of v (here we dropped the some- 
what arbitrary factor -y/Tr in (|50p ).The weakly nonlinear 
approximation holds to the south-west of these curves 
when ^/ew <C 1. The criterion given by the inequality 
([5T|) corresponds to the asymptotic behaviour of those 
curves as v tends to infinity. 

To observe that this is the case numerically, we have 
also compared the charging dynamics given by the weakly 
nonlinear approximation to that of the full MPNP solu- 
tion in Figini with the parameters e — v = 0.01. When 
the parameters were to the south-west of the = 1, 
as in case (c), the match was perfect. In the other cases, 
the weakly nonlinear approximation did not do as well, it 
was particularly off for case (b), when the product ^/Iw 
yielded the highest number. In case (a), the agreement 
was off by a constant shift, whereas in case (d), the situa- 
tion improved over time. This is because ew is still small 
for case (d), although ^/ew is not, and therefore weakly 
nonlinear approximation is still valid for the final state 
of the system. To summarize, if an accurate description 
of the system at all times is desired, then ^/ew << 1 is 
the appropriate criterion, however it may suffice to have 
just ew << 1 to be able to predict the eventual steady 
state by the weakly nonlinear model. 



V. CONCLUSIONS 

As an extension of the modified PB approach, we have 
derived modified PNP equations, which may be of help 
when the thin double layer approximation fails. Using 
this new set of equations, we have revisited the model 
problem of the step charging of a cell with parallel block- 
ing electrodes. In addition, we have confirmed through 
asymptotic analysis the hypotheses stated in Part I re- 
garding the MPB double layer model. We have also inves- 
tigated the limits of the thin double layer approximation 
as well as higher order corrections. 

The MPNP system proposed is a natural extension 
of one of the two models introduced in Part I, namely 
the MPB model based on an approach originally due 
to [i^, HSl. One can similarly construct other MPNP 
equations from other MPB models such as the composite 



diffuse-layer model also introduced in Part I. However the 
discontinuous structure of the latter leads to a complex 
formalism with discontinuities, improper for implemen- 



Curves of e^'^w = 1 




20 40 60 80 100 

M (zeV/kT) 



FIG. 4: Shown are the curves on which ^/ew — 1 for various 
values of v. The weakly nonhnear approximation holds to the 
south-west of these curves when ^/ew << 1. 



tation in complex geometries, whereas the one presented 
here has a smooth behavior and is therefore much more 
broadly applicable. 

We expect that the MPNP equations presented here, 
or other simple variations with different modifications of 
the chemical potentials or free energy, will find many 
applications. They are no more difficult to use than the 
classical PNP equations, which are currently ubiquitous 
in the modeling of electrochemical systems. Especially 
at large voltages, the MPNP equations are much better 
suited for numerical computations, since they lack the 
exponentially diverging concentrations predicted by the 
PNP equations, which are difficult to resolve. Of course, 
those same divergences are also clearly unphysical, while 
the MPNP equations predict reasonable ion profile for 
any applied voltage. 

Future research directions include the application of 
the presented framework to many other settings including 
the response of the simple electrolytic cell considered here 
to various systems with time-dependent applied voltages 
of strong amplitudes. In particular, driving an electro- 
chemical cell or AC electro-osmotic pump at rather large 
frequencies should be a selective way of checking the va- 
lidity/use of equations such as the one put forward here 
because dynamical effects are exacerbated in such situa- 
tions. 
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FIG. 5; Four case studies on the validity of the weakly nonlin- 
ear approximation. No stern capacitance is included ((5 = 0). 
(I) The y/ew = 1 curve for the particular value f — 0.01, and 
corresponding locations of the four case studies shown. (II) 
The comparison of the charge q{t) stored in the half cell com- 
puted by (i) the MPNP (ii) the weakly nonlinear approxima- 
tion corresponding to the MPNP. The match in (c) is good as 
expected, and the match in (b) is off by several factors, again 
as expected. In cases (a) and (d), the curves run close but 
they are clearly separated. 
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